Determining the stability of genetic switches: exphcitly accounting for mRNA noise 
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Cells use genetic switches to shift between alternate gene expression states, e.g., to adapt to new 
environments or to follow a developmental pathway. Here, we study the dynamics of switching in a 
generic-feedback on/off switch. Unlike protein-only models, we explicitly account for stochastic fluc- 
tuations of mRNA, which have a dramatic impact on switch dynamics. Employing the WKB theory 
to treat the underlying chemical master equations, we obtain accurate results for the quasi-stationary 
distributions of mRNA and protein copy numbers and for the mean switching time, starting from 
either state. Our analytical results agree well with Monte Carlo simulations. Importantly, one can 
use the approach to study the effect of varying biological parameters on switch stability. 
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Genetic switches allow cells to switch between distinct 
gene expression states in response to environmental stim- 
uli and/or internal signals. The ultimate stability of 
these states is determined by stochastic fluctuations of 
mRNA and proteins during gene expression [T] that can 
give rise to spontaneous switching, even in the absence of 
a driving signal. When gene expression states are stable 
on the time scale of cellular division they can carry epige- 
netic information across generations, however, when they 
are more transient they may provide a beneficial source 
of heterogeneity in genetically identical populations. 

Previous studies of noise-driven genetic switches have 
shown that switching can be treated as a first-passage 
problem of the underlying Markov process. These results 
were obtained either by using the probability generating 
function formalism, or by employing semi-classical ap- 
proximation schemes to the chemical master equations 
(CMEs) or related Langevin equations [IHSj. These stud- 
ies, however, focused on protein-only models and ignored 
the presence of mRNA and thereby the influence of tran- 
scriptional noise. Recently it has been shown that ex- 
plicitly accounting for mRNA in the underlying CMEs 
has a strong impact on switching times O [7] . A general 
method for accurately determining the mRNA/protein 
distributions and stability of feedback-based switches is 
of great interest, as they regulate diverse biological phe- 
nomena, such as microbial environmental adaptation, de- 
velopmental pathways, and bacteriophage lysogeny [THS]. 

In this study, we explicitly account for the mRNA noise 
and present a concise analytical framework for accurately 
calculating the stability of gene expression switches sub- 
ject to stochastic fluctuations of protein/mRNA. Our ap- 
proach is demonstrated on a two-state positive feedback 
switch, which was experimentally shown to describe bi- 
ological switching [TUl. We apply a WKB theory [TTl 
to the CMEs and obtain the quasistationary probabil- 
ity distributions of mRNA and protein copy numbers in 
the on and off states, from which we extract the mean 
switching times starting from either state. Our results 
agree well with Monte Carlo simulations. Finally, we use 
our analytical predictions to study the effect of promoter 
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FIG. 1: (Color online) (a) Model for positive feedback net- 
work. Transcription and translation are modeled as first-order 
processes with rates a and 76, respectively. mRNA and pro- 
teins undergo first order degradation with rates 7 and 1 (we 
rescale rates by the protein decay rate). The feedback func- 
tions kon{n) and koff{n) control promoter transitions, (b) 
The momentum py vs protein copy number n. The thick line 
indicates the off^on switching trajectory, and shaded areas 
correspond to the entropy barriers for switching, (c) mRNA 
m and protein counts in a typical Monte Carlo trajectory un- 
dergoing switching. In (b,c) K^ab = 2400, b = 22.5, h^2, 
nso^lOOO, fcr" = fer'" = a/100, fcj^"" = fcr" =«, and 7 = 50. 



fluctuations on switching stability. 

We consider a two-state gene-expression model where 
transitions between a transcriptionally active and inac- 
tive promoter are controlled by the protein copy number 
n via positive feedback (see Fig. [TJa)). The transition 
rates into the active and inactive states are ko„in) = 
f{n) and fco/^(n) = f/(n). While our analytical treatment 
holds for generic f{n) and g{n), we consider a con- 
crete example using Hill- type functions /(n) = 
(/cJJ^'»^-fc^")n''V(n5(5 +n^^) and g{n) = kY''''' -{k'P'''' - 
fc™™)n''^/(n5Q -I- n^^), which were shown to be biologi- 
cally relevant, e.g., in the lac operon j^. Here nso is the 
curve's midpoint and for simplicity we set hi = h2 = h. 
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FIG. 2: (Color online) (a) Protein (left) and mRNA (right) 
QSDs in the off state showing WKB result ^ (solid) and 
MC simulations (A) for b — 2 and h — 2. Our results converge 
to those of |12) (dashed) for h ^ ex. (b) As in (a) for the 
on state, (c) The KuUback-Leibler divergence (see text) vs. 
h comparing the WKB and MC PDFs of the off state for 
6 = 2.5. (d) The KL divergence vs. b of the off state for h — 2 
(x) and /i = 2.5 (A). Insets: The theory holds for KAS > 1. 
Other parameters are 7^ = 06 = 3200, n5o = 2000, k^'" = a/50, 
A;n"=a, kr" = a/100, fcr''" = a/2, and 7 = 50. 

The deterministic rate equations (DREs) for the mean 
number of mRNA, M, and proteins, N, read 

M = afiN)/[f{N)+g{N)]~jM, N^jbM-N, (1) 

where f{N)/[f{N) + g{N)] is the probability for an ac- 
tive promoter. To exhibit bistabihty, Eqs. ([ij must have 
(at least) three (positive) fixed points. We denote by 
Non and Noff, respectively, the attracting fixed points 
corresponding to the average protein copy number in the 
on and off states, and assume that 1 <^ ^off ^ Non- 
These points are separated by a repeller A^o such that 
Noff < Nq < Non- For biologically-relevant parameters, 
see below, one has a - /c^f^ > and 5 = 0(1). 

Also, when h 1, Non ^ ab NoffkY''''=/kJJ''". Thus, 
K = ab^ 1 - the typical protein number in the on state 
- will serve as the large parameter of the theory. 

DREs ^ ignore noise and predict that, once the sys- 
tem has settled in one of the attracting fixed points, it 
stays there forever. Yet, the presence of intrinsic noise al- 
lows switching between these fixed points by crossing the 
corresponding entropy barrier [T^] . In the stochastic pic- 
ture, starting from the vicinity of either state the system 
rapidly converges into the quasistationary distribution 
(QSD) about this state. This distribution is metastable, 
and slowly decays due to a (exponentially small) probabil- 
ity leakage through the entropy barrier at A^o O El] ■ It 
is this leakage that determines the corresponding switch- 
ing rates between the metastable states. 

To model the stochastic behavior of the switch, we use 



two coupled CMEs. These describe the dynamics of Pm.n 
and Qm,n ~ the probability distribution functions (PDFs) 
of having m mRNAs and n proteins at time t with the 
promoter in the inactive and active state, respectively: 

Qm,n = -9{n)Qm,n + I{n)P„,,n + [A+a(£;-i - 1)] Q™,„.(2) 

Here, EU{n)^^ f{n + j), A = {Ei - l)n + ^{E^ - 
l)m + 'ybm{Ej^ ^ — I) is a birth-death operator related to 
the inactive promoter, and ^ Pm,n + Qm,n — 1- We are 
seeking the QSD starting from the vicinity of Noff (the 
on state treatment is equivalent, see below). Putting 
Pm,n = Qm.n = in Eqs. ^ (that are exponentially 
small for K ^ 1), and eliminating, e.g., Qm,n we obtain 

0={A + 5(n)-i [A + a(£;„i-I)][/(n)-A]}P™,„. (3) 

In the case of constant transition rates between the 
active and inactive states, Eqs. ([3| were asymptotically 
solved in the 7^1 limit using the probability generating 
function [T^[TS]. However, for generic feedback functions 
the generating function formalism cannot be used. 

Instead, we use here a powerful method based on the 
WKB approximation [TT], to treat the (quasi)stationary 
CMEs ^ [13^. The WKB ansatz reads 

P,n,n = P{x, y) ^ eM-KS{x, y)] . (4) 

Here x = m/K and y = n/K are the densities of the 
mRNA and proteins, respectively. Plugging ansatz ^ 
into Eq. ([s]), e.g. the step operator E^ is replaced in 
the leading order by the function ef'^^^i^-v) , After some 
algebra, one arrives in the leading order at a stationary 
Hamilton-Jacobi equation H(x,y,dxS,dyS) — 0, with 

H^A+~giyr' [A + b'\eP- - 1)] [f{y)-A]. (5) 

Here A = A{x,y,Px,Py) = y{e~Py-l) + 7a;(e"P- -1) -f- 
"fbx{ePy — 1) is now a function, and we have used the 
rescaled feedback functions f{y) — f{y)/K and g{y) = 
g{y)/K, and mRNA production rate a/K = b~^. Also, 
in analogy to classical mechanics we have introduced the 
momenta p^ = dxS{x,y) and Py = dyS{x,y), corre- 
sponding to the steepness of the sought PDF's [T?l I16j . 
Note that putting p^ ~ py = corresponds to mean-field 
dynamics; in this case using ^ the Hamilton's equations 
for X = dp^ H and y = dp^ H become DREs ([T]) divided by 
.9('^)/[/('^) + .9('^)] [since ^ corresponds to the inactive 
state occupied with probability g{n)/[f{n) + g{n)]]. 

The strength of this theory is that it can accurately ac- 
count for rare large fluctuations responsible for switching. 
To do so one has to solve the Hamilton's equations for x 
and y together with p^ — —dxH and Py — —dyH. Now, 
as we look for a zero-energy trajectory of ([5]), the action 
reads S{x,y) = J p^dx+pydy, which yields PDF 2D 
Hamiltonian systems can be solved numerically jS^. In 



3 



our case such a solution would yield the complete statis- 
tics for arbitrary mRNA degradation rates, which is im- 
portant e.g. in eukaryotic systems where 7 = C(l). 

Further analytical progress can be made in the regime 
of 7 ^ 1 (relevant for bacterial systems), for which the 
mRNA dynamics is enslaved to that of the protein [T^. 
We adiabatically eliminate [T7] the fast component in 
the mRNA dynamics by assuming x and Px rapidly con- 
verge to slowly- varying functions of {y,Py). Taking y 
and Py constant and putting x = p^ = 0, one obtains 
X = 0(7-1) ^8j, andp:, = -ln{l + b-bePy ) [H], so that 
Eq. (l5| becomes a reduced Hamiltonian Hr{y,Py) 
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This Hamiltonian effectively accounts for 



the fact that the proteins are produced in geometrically 
distributed bursts with mean &, which in turn asymptot- 
ically accounts for the mRNA noise when 7 1. How- 
ever, had one initially eliminated the mRNA species from 
CMEs ([2]) using geometrically distributed protein births, 
one would have obtained an analytically intractable one- 
dimensional CME. This is because the system under con- 
sideration is a two-state switch with non-linear feedback. 

The (nontrivial) zero-energy trajectory of Hamilto- 
nian ([6| encodes the stochastic dynamics of (only) the 
proteins, and corresponds to its quasi-stationary behav- 
ior. This trajectory gives Py as function of y and repre- 
sents the most probable path the stochastic system fol- 
lows while undergoing switching |13l 116) . The normaliz- 
able solution reads pj,(y) = \n[{-B+y/W^^lAC) / {2A% 
see Fig. 0b), where A ^ {l + by)[y + j{y)]+byg[y), B = 
-y[yil + 2b) + l + {l + b){f{y) + ~g{y))], and C = {l + b)y^. 
Thus S{y) ~ Pviy')dy' |20j, and using Q we have 
P{y) ^ e~^^'^y\ Note that P{y) is the contribution to 
the QSD corresponding to an inactive promoter. A sim- 
ilar contribution from the active prornoter can be shown 
to also satisfy Q{y) ~ e 



tein copy number QSD, Vn, starting from the vicinity of 
the off state, reads Vn = V{y) = P{y)+Q{y) ~ e'-^^^^'). 
Expanding S{y) in the vicinity of yo/f — No/f/K up to 
second order, and demanding that the Gaussian integra- 
tion be normalized to 1, the normalized V{y) satisfies 



V{y) ^ y^5"(y„;/)/(27ri^)e-^[^(^)-^(''°//)l. (7) 

Note that the preexponent entering ([?]) holds only in 
the Gaussian regime of the PDF, whose width is cr = 
^/K/S"{yoff). One can check that the on state QSD 
coincides with ([7| upon replacing yoff~^yon- The Gaus- 
sian normalization above is valid when the QSD's width 
is sufficiently small compared to Noff for the off state 
(and Non — Nq for the on state) . From ([7| one can read- 
ily find the joint QSD, Vm,n = 'Pm\n'Pn, where Vm\n, the 
probability to find m mRNA molecules given n proteins. 
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FIG. 3: (Color online) (a) MST Tojf^on as a function of b for 
nso — 720. WKB result ([8| with numerical preexponent (solid) 
and MC simulations (xjT^(b) Ton^off as a function of nso for 
b = 15. Other parameters are h = 2, K = at = 2400, /fc™"' = 
fcr"' = a/100, fcj'''^ = fer™=a and 7 = 50. Preexponents were 
17.8 (a) and 4.8 (b). Also shown are the MSTs from MC 
simulations of protein burst models with constant (Q) and 
geometrically distributed (□) burst sizes. Insets: WKB result 
without (dashed) and with (solid) numerical preexponent. 



can be found using standard techniques |12| . Given Vm.m 
the mRNA QSD satisfies Vm = E„ 'Pm,n- 

Our result ([7| can be compared to that of Ref. [12] e.g. 
for the one-state model, where the promoter is always 
active. Putting f{y) = 1 and g{y) = 0, the momentum 
becomes py{y) = ln[{{l+b)y)/{l+by))] so that Eq. Q sim- 
plifies to Vn = [2TTab{b+l)]-^^'^a-''n-"'{a + n)''+'''V'{l + 
This result coincides with the n ^ 1 asymp- 
tote of Eq. (9) in [T^] by using the Stirling formula. As 
expected, the prefactor here coincides with that of [12] 
only in the Gaussian region of the fixed point n = ab. 

To check our theoretical predictions for generic non- 
constant f{n) and g{n), we performed Monte Carlo (MC) 
simulations using the Gillespie algorithm [21 . An exam- 
ple of a typical MC run can be seen in Fig. [ijc). In 
Fig. [2][a,b) we compare the WKB prediction (where the 
action is found by numerical integration) for the protein 
and mRNA QSDs for the off (a) and the on (b) states, 
with MC simulations and results of Ref. [T^. The lat- 
ter are expected to be valid only in the limit of /i ^ 1 
(when the feedback functions become approximately step 
functions). In panels (c,d) we show the Kullback-Leibler 
(KL) divergence ^ P^^^ ln(Pi"^^/Pi^^) (a measure of the 
difference between PDFs Pn^^ and P,!^'') between WKB 
result ([7| and MC simulations for various parameters. 

Now, the mean switching time (MST) Toff^on can be 
readily inferred from QSD Q: it is the inverse of the 
flux through the repelling flxed point yQ^N^/K [T31 [11] . 
The logarithm of the MST is proportional to the effective 
entropy barrier between the attracting and repelling fixed 
points AS'o// = 5'(yo) — S{yoff). Therefore, we have 
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= K[ASoff + 0{l/K,l/j)] 
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with T, 



on^off 



the same using AS'o 



S{yo) - S{yo 



Note, that while the prefactor of the MST is unknown, 
based on single-species calculations it is expected to be 
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ized by metastable switching, e.g., those with additional 
promoter states such as DNA looping or nucleosome re- 
modeling. In particular, it can be used to help elucidate 
the underlying regulatory circuits responsible for pheno- 
typical changes as a result of switching. 
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FIG. 4: (Color online) (a) To. 



off /Tof f^on (on state relative 
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fcr"' = aQ, if = 2400, 6=15, andn5o = 950. (b) MSTs r<,„_+o// 
(dashed), Toff^on (dash-dotted), and their ratio (solid). 



stability) vs h and a, for fcj'™ = k^" = aa/100, 



0(1) and independent of K [TB]. Eq. ^ indicates the 
WKB formalism is valid for KAS ^ 1 (see Fig. [2] insets). 

In Fig. [3] we compare the theoretical MST prediction 
to MC simulations. Panel (a) compares Toff^on vs b; a 
nontrivial super-exponential dependence is observed. It 
is also shown that the functional dependence of Tq/ f^on 
is excellently captured by ([s]) with a numerical slowly- 
varying preexponent. Panel (b) compares Ton^off vs nso, 
and again the functional dependence is excellently cap- 
tured by the theoretical MST. Fig. [3] also shows that the 
mRNA noise can be accounted for when 7 ^ 1 by assum- 
ing geometrically distributed protein births, while using 
a constant burst size yields markedly different results. 

Promoter fluctuations are known to have a substan- 
tial impact on the PDFs and switching times of genetic 
switches [H |5] . We used Eq. ([s]) to study switch stability 
with respect to the promoter transition dynamics. Mul- 
tiplying kon and fco// by a we control the frequency and 
duration of mRNA bursts in the off state and pauses in 
the on state, while leaving unchanged the relative proba- 
bility of the promoter to be in the active/inactive states. 

By increasing a the duration of bursts and pauses di- 
minishes, large fluctuations become rarified, and crossing 
the entropic barrier becomes harder (see also [S]). The 
rate of increase of the entropic barrier, however, differs 
for the off — ?> on and on — ^ off transitions, as seen in 
Fig. |4] For small h the rate of increase of the on — > off 
entropic barrier exceeds that of off — > on for the en- 
tire range of a, thereby amplifying the on state stability. 
However, this is not the case for higher h, which gives 
rise to a non-monotonic stability curve for the on state 
(see panels b in Fig. |4| . These results stress the role of 
promoter kinetics, not just thermodynamics, for genetic 
switches, which are inherently far from equilibrium. 

We have presented an analytical framework for the ac- 
curate analysis of genetic switches while explicitly ac- 
counting for mRNA noise. This framework is expected to 
be useful for studying diverse genetic circuits character- 
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